* Figure and table: Bayesian estimation

use "${output}panel_r0_r1_r2.dta", clear // Load merged analysis panel

********************************************************************************

* Keep households surveyed in each round
keep if sample_household == 1

********************************************************************************

* Unique sort 
sort hh_id surveyround

* Rename ICS purchase variable
ren purchased_intervention_stove y

* MCMC specification
local burn_in = 10000
local mcmcsize = 50000
local thinning = 5

*********************************************************************************

* INFORMATIVE PRIOR

local mu = -1.17
local sigma2 = 1

// Generate district dummies
tab districtcode, gen(districtdummy)
ren districtcode discode

reghdfe y i.treatment##i.chirag_strata if surveyround == 1, absorb(discode) vce(cluster uniquegrp)
local treatment_prior_mean = _b[1.treatment]
local chirag_prior_mean = _b[1.chirag_strata]
local treatment_chirag_prior_mean = _b[1.treatment#1.chirag_strata]

bayes, prior({y:1.treatment#1.chirag_strata}, lognormal(`mu',`sigma2')) ///
	prior({y:1.treatment}, normal(`treatment_prior_mean',10000)) ///
	prior({y:1.chirag_strata}, normal(`chirag_prior_mean',10000)) ///
	rseed(17642) ///
	thinning(`thinning') ///
	mcmcsize(`mcmcsize') burnin(`burn_in') ///
	saving("${output}inform_simulation_results.dta", replace): ///
		mixed y i.treatment##i.chirag_strata if surveyround == 1 || discode: || gp_code: || uniquegrp:

// Extract parameters for Bayes table
foreach x in 1.treatment 1.chirag_strata 1.treatment#1.chirag_strata {

	// Adjust name for local macros
	local y `x'
	local y = subinstr("`y'",".","",.)
	local y = subinstr("`y'","#","",.)
	
	// Extract from table
	local inf_m_`y' = e(mean)[1,colnumb(e(mean),"`x'")] // mean
	local inf_sd_`y' = e(sd)[1,colnumb(e(sd),"`x'")] // sd
	local inf_l_`y' = e(cri)[1,colnumb(e(cri),"`x'")] // lower credible interval
	local inf_u_`y' = e(cri)[2,colnumb(e(cri),"`x'")] // upper credible interval

	foreach i in m sd l u {
		local inf_`i'_`y' : di %5.2f `inf_`i'_`y'' // 2 decimal places

	}
}

bayesstats ess

bayesgraph diagnostics {y:1.treatment#1.chirag_strata}
graph export "${output}inform_bayes_diagnostics_interaction.png", replace width(4000)

bayesgraph matrix {y:1.treatment#1.chirag_strata} {y:1.treatment} {y:1.chirag_strata}
graph export "${output}inform_bayes_correlation.png", replace width(4000)

* Graph of prior and posterior distribution

preserve

	use "${output}inform_simulation_results.dta", clear

	sum eq7_p8, detail // Check that this is the saved MCMC sample of estimates for TREATMENT X NGO
	
	graph twoway (kdensity eq7_p8, 	///
			xscale(range(0 1)) ///
			yscale(range(0 6)) ylabel(0(1)6)) ///
		|| (kdensity eq7_p8, 	///
			range(`r(p5)' `r(p95)') color(gs14%90) recast(area) lwidth(0)) ///
		|| (function y = (1 / x) * normalden(ln(x),`mu',`sigma2'), range(0 1) ///
			ytitle("") xtitle("") ///
			title("") ///
			lpattern(dash) ///
			title("(a) Main result: Informative prior") ///
			legend(order(3 1 2) label(3 "Prior") label(1 "Posterior") label(2 "90% credible interval") ///
			ring(0) pos(1)))
		
	graph export "${output}inform_bayes_prior_posterior.png", replace width(4000)
	graph save "${output}inform_bayes_prior_posterior.gph", replace

restore

*********************************************************************************

* NON-INFORMATIVE PRIOR

bayes, prior({y:1.treatment}, normal(`treatment_prior_mean',10000)) ///
	prior({y:1.chirag_strata}, normal(`chirag_prior_mean',10000)) ///
	prior({y:1.treatment#1.chirag_strata}, normal(`treatment_chirag_prior_mean',10000)) ///
	rseed(17642) ///
	thinning(`thinning') ///
	mcmcsize(`mcmcsize') burnin(`burn_in') ///
	saving("${output}noninform_simulation_results.dta", replace): ///
		mixed y i.treatment##i.chirag_strata if surveyround == 1 || discode: || gp_code: || uniquegrp:

// Extract parameters for Bayes table
foreach x in 1.treatment 1.chirag_strata 1.treatment#1.chirag_strata {

	// Adjust name for local macros
	local y `x'
	local y = subinstr("`y'",".","",.)
	local y = subinstr("`y'","#","",.)
	
	// Extract from table
	local dif_m_`y' = e(mean)[1,colnumb(e(mean),"`x'")] // mean
	local dif_sd_`y' = e(sd)[1,colnumb(e(sd),"`x'")] // sd
	local dif_l_`y' = e(cri)[1,colnumb(e(cri),"`x'")] // lower credible interval
	local dif_u_`y' = e(cri)[2,colnumb(e(cri),"`x'")] // upper credible interval

	foreach i in m sd l u {
		local dif_`i'_`y' : di %5.2f `dif_`i'_`y'' // 2 decimal places
		
	}
}

bayesstats ess

bayesgraph diagnostics {y:1.treatment#1.chirag_strata}
graph export "${output}noninform_bayes_diagnostics_interaction.png", replace width(4000)

bayesgraph matrix {y:1.treatment#1.chirag_strata} {y:1.treatment} {y:1.chirag_strata}
graph export "${output}noninform_bayes_correlation.png", replace width(4000)

* Graph of prior and posterior distribution

preserve

	use "${output}noninform_simulation_results.dta", clear

	sum eq7_p8, detail // Check that this is the saved MCMC sample of estimates for TREATMENT X NGO

	graph twoway (kdensity eq7_p8) 	///
		|| (kdensity eq7_p8, 	///
			range(`r(p5)' `r(p95)') color(gs14%90) recast(area) lwidth(0)) ///
		|| (function y = normalden(x,0,10000), range(-0.2 0.5) ///
		ytitle("") xtitle("") ///
		lpattern(dash) ///
		title(`"(b) Robustness: Diffuse ("uninformative") prior"') ///
		legend(order(3 1 2) label(3 "Prior") label(1 "Posterior") label(2 "90% credible interval") ///
		ring(0) pos(1)))
		
	graph export "${output}noninform_bayes_prior_posterior.png", replace width(4000)
	graph save "${output}noninform_bayes_prior_posterior.gph", replace

restore

*********************************************************************************

* STRONG PRIOR ON NGO-REPUTATION-EFFECT = 0

bayes, prior({y:1.treatment#1.chirag_strata}, normal(0, 0.01)) ///
	prior({y:1.treatment}, normal(`treatment_prior_mean',10000)) ///
	prior({y:1.chirag_strata}, normal(`chirag_prior_mean',10000)) ///
	rseed(17642) ///
	thinning(`thinning') ///
	mcmcsize(`mcmcsize') burnin(`burn_in') ///
	saving("${output}strongnull_simulation_results.dta", replace): ///
		mixed y i.treatment##i.chirag_strata if surveyround == 1 || discode: || gp_code: || uniquegrp:
		
// Extract parameters for Bayes table
foreach x in 1.treatment 1.chirag_strata 1.treatment#1.chirag_strata {

	// Adjust name for local macros
	local y `x'
	local y = subinstr("`y'",".","",.)
	local y = subinstr("`y'","#","",.)
	
	// Extract from table
	local str_m_`y' = e(mean)[1,colnumb(e(mean),"`x'")] // mean
	local str_sd_`y' = e(sd)[1,colnumb(e(sd),"`x'")] // sd
	local str_l_`y' = e(cri)[1,colnumb(e(cri),"`x'")] // lower credible interval
	local str_u_`y' = e(cri)[2,colnumb(e(cri),"`x'")] // upper credible interval

	foreach i in m sd l u {
		local str_`i'_`y' : di %5.2f `str_`i'_`y'' // 2 decimal places
		
	}
}

bayesstats ess

bayesgraph diagnostics {y:1.treatment#1.chirag_strata}
graph export "${output}strongnull_bayes_diagnostics_interaction.png", replace width(4000)

bayesgraph matrix {y:1.treatment#1.chirag_strata} {y:1.treatment} {y:1.chirag_strata}
graph export "${output}strongnull_bayes_correlation.png", replace width(4000)

preserve

	use "${output}strongnull_simulation_results.dta", clear

	sum eq7_p8, detail // Check that this is the saved MCMC sample of estimates for TREATMENT X NGO

	graph twoway (kdensity eq7_p8) 	///
		|| (kdensity eq7_p8, 	///
			range(`r(p5)' `r(p95)') color(gs14%90) recast(area) lwidth(0)) ///
		|| (function y = normalden(x,0,0.01), range(-0.2 0.6) ///
		ytitle("") xtitle("") ///
		lpattern(dash) ///
		title("(c) Robustness: Strong no-NGO-reputation-effect prior") ///
		legend(order(3 1 2) label(3 "Prior") label(1 "Posterior") label(2 "90% credible interval") ///
		cols(3) pos(6)))
		
	graph export "${output}strongnull_bayes_prior_posterior.png", replace width(4000)
	graph save "${output}strongnull_bayes_prior_posterior.gph", replace

restore

*********************************************************************************

* Create graphs for manuscript 

gr drop _all
grc1leg ///
	"${output}inform_bayes_prior_posterior.gph" ///
	"${output}noninform_bayes_prior_posterior.gph" ///
	"${output}strongnull_bayes_prior_posterior.gph", ///
	legendfrom("${output}strongnull_bayes_prior_posterior.gph") ///
	b1("{&beta}(TREATMENT{sub:j} × NGO{sub:j})") l1("Density") ///
	name(bayesian_combined, replace) cols(1) ring(1) position(6)

gr display bayesian_combined, ysize(10)

graph export "${results}figure_bayesian_estimation_main_robustness.png", width(5000) replace
gr drop _all

*********************************************************************************

* Generate appendix results table

texdoc init "${results}apptable_bayesian_estimation.tex", replace force
	tex \begin{table}[tbph]
	tex \centering
	tex \caption{Markov Chain Monte Carlo results}
	tex \begin{adjustbox}{max width=\textwidth}
	tex \label{app:tab_mcmc}
	tex \begin{threeparttable}
	tex \begin{tabular}{lccc} \toprule 
	tex & (1) & (2) & (3) \\
	tex & Posterior mean & Posterior standard deviation & 95\% credible interval \\ \midrule \addlinespace
	tex \multicolumn{4}{l}{\textbf{(a) Model I}} \\
	tex \addlinespace
	tex \({TREATMENT}_j\) & `inf_m_1treatment' & `inf_sd_1treatment' & \(\left[`inf_l_1treatment', `inf_u_1treatment' \right]\) \\
	tex \({NGO}_j\) & `inf_m_1chirag_strata' & `inf_sd_1chirag_strata' & \(\left[`inf_l_1chirag_strata', `inf_u_1chirag_strata'\right]\) \\
	tex \({TREATMENT}_j \times {NGO}_j\) & `inf_m_1treatment1chirag_strata' & `inf_sd_1treatment1chirag_strata' & \(\left[`inf_l_1treatment1chirag_strata', `inf_u_1treatment1chirag_strata'\right]\) \\
	tex \addlinespace
	tex \midrule \addlinespace
	tex \multicolumn{4}{l}{\textbf{(b) Model II}} \\
	tex \addlinespace
	tex \({TREATMENT}_j\) & `dif_m_1treatment' & `dif_sd_1treatment' & \(\left[`dif_l_1treatment', `dif_u_1treatment' \right]\) \\
	tex \({NGO}_j\) & `dif_m_1chirag_strata' & `dif_sd_1chirag_strata' & \(\left[`dif_l_1chirag_strata', `dif_u_1chirag_strata'\right]\) \\
	tex \({TREATMENT}_j \times {NGO}_j\) & `dif_m_1treatment1chirag_strata' & `dif_sd_1treatment1chirag_strata' & \(\left[`dif_l_1treatment1chirag_strata', `dif_u_1treatment1chirag_strata'\right]\) \\
	tex \addlinespace
	tex \midrule \addlinespace
	tex \multicolumn{4}{l}{\textbf{(c) Model III}} \\
	tex \addlinespace
	tex \({TREATMENT}_j\) & `str_m_1treatment' & `str_sd_1treatment' & \(\left[`str_l_1treatment', `str_u_1treatment' \right]\) \\
	tex \({NGO}_j\) & `str_m_1chirag_strata' & `str_sd_1chirag_strata' & \(\left[`str_l_1chirag_strata', `str_u_1chirag_strata'\right]\) \\
	tex \({TREATMENT}_j \times {NGO}_j\) & `str_m_1treatment1chirag_strata' & `str_sd_1treatment1chirag_strata' & \(\left[`str_l_1treatment1chirag_strata', `str_u_1treatment1chirag_strata'\right]\) \\
	tex \bottomrule
	tex \end{tabular}
	tex \begin{tablenotes}
	tex {\setlength\labelsep{0pt}
	tex \footnotesize
	tex \item \textit{Notes}. The outcome variable is an indicator that equals 1 if household \(i\) in hamlet \(j\) purchased at least one of the two ICS promoted during the intervention. Columns (1) and (2) present the mean and standard deviations, respectively, for the MCMC sample. Column (3) presents the 95 percent credible interval for the MCMC sample. Estimates for the district fixed-effect and the 97 hamlet-specific random effects not reported for brevity. As in Table \ref{t:purchase_rates}, \(N = 943\) households.}
	tex \end{tablenotes}
	tex \end{threeparttable}
	tex \end{adjustbox}
	tex \end{table}
texdoc close

*********************************************************************************
